Fit temperature models and predict growing season temperature

Author

Max Lindmark, Jan Ohlberger, Anna Gårdmark

Published

July 6, 2023

Load libraries

pkgs <- c("here","tidyverse", "tidylog", "RColorBrewer", "viridis", "sdmTMB", "sdmTMBextra", "patchwork", "RCurl", "tidylog") 

# minpack.lm needed if using nlsLM()
if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
  }

invisible(lapply(pkgs, library,character.only = T))
here() starts at /Users/maxlindmark/Dropbox/Max work/R/perch-growth
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.2     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.2.1
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Attaching package: 'tidylog'


The following objects are masked from 'package:dplyr':

    add_count, add_tally, anti_join, count, distinct, distinct_all,
    distinct_at, distinct_if, filter, filter_all, filter_at, filter_if,
    full_join, group_by, group_by_all, group_by_at, group_by_if,
    inner_join, left_join, mutate, mutate_all, mutate_at, mutate_if,
    relocate, rename, rename_all, rename_at, rename_if, rename_with,
    right_join, sample_frac, sample_n, select, select_all, select_at,
    select_if, semi_join, slice, slice_head, slice_max, slice_min,
    slice_sample, slice_tail, summarise, summarise_all, summarise_at,
    summarise_if, summarize, summarize_all, summarize_at, summarize_if,
    tally, top_frac, top_n, transmute, transmute_all, transmute_at,
    transmute_if, ungroup


The following objects are masked from 'package:tidyr':

    drop_na, fill, gather, pivot_longer, pivot_wider, replace_na,
    spread, uncount


The following object is masked from 'package:stats':

    filter


Loading required package: viridisLite


Attaching package: 'sdmTMBextra'


The following objects are masked from 'package:sdmTMB':

    dharma_residuals, extract_mcmc



Attaching package: 'RCurl'


The following object is masked from 'package:tidyr':

    complete
# devtools::install_github("seananderson/ggsidekick") # not on CRAN 
library(ggsidekick)
theme_set(theme_sleek())

# Set path:
home <- here::here()

Load cache

# qwraps2::lazyload_cache_dir(path = paste0(home, "/R/analyze-data/01-fit-temp-models-predict_cache/html"))

Read data

d <- read_csv(paste0(home, "/output/temp_data_for_fitting.csv"))
Rows: 98616 Columns: 13
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (4): area, source, db, id
dbl  (8): year, temp, yday, month, day, VkDag, stn_nr, section_nr
date (1): date

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
d <- d |> mutate(area = as.factor(area),
                 source_f = as.factor(source),
                 year_f = as.factor(year))
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
# Keep track of the years for which we have cohorts for matching with cohort data
gdat <- readr::read_csv("https://raw.githubusercontent.com/maxlindmark/perch-growth/master/data/for-analysis/dat.csv")
Rows: 338460 Columns: 12
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (6): age_ring, area, gear, ID, sex, keep
dbl (6): length_mm, age_bc, age_catch, catch_year, cohort, final_length

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
gdat$area_cohort_age <- as.factor(paste(gdat$area, gdat$cohort, gdat$age_bc))

gdat <- gdat |>
  group_by(area_cohort_age) |> 
  filter(n() > 10) |> 
  filter(age_catch > 3) |> 
  group_by(area) |>
  summarise(min = min(cohort),
            max = max(cohort)) |> 
  arrange(min)
group_by: one grouping variable (area_cohort_age)
filter (grouped): removed 5,190 rows (2%), 333,270 rows remaining
filter (grouped): removed 107,058 rows (32%), 226,212 rows remaining
group_by: one grouping variable (area)
summarise: now 12 rows and 3 columns, ungrouped
d <- left_join(d, gdat, by = "area") |>
  mutate(area = as.factor(area),
         growth_dat = ifelse(year >= min & year <= max, "Y", "N"))
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (     0)
           > matched rows     98,616
           >                 ========
           > rows total       98,616
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
# Drop data in SI_HA and BT before onset of warming
d <- d |>
  mutate(discard = "N",
         discard = ifelse(area == "BT" & year <= 1980, "Y", discard),
         discard = ifelse(area == "SI_HA" & year <= 1972, "Y", discard)) |> 
  filter(discard == "N")
mutate: new variable 'discard' (character) with 2 unique values and 0% NA
filter: removed 1,146 rows (1%), 97,470 rows remaining
# Drop heated areas actually for the full plot
df <- d |> filter(!area %in% c("BT", "SI_HA"))
filter: removed 24,149 rows (25%), 73,321 rows remaining

Fit models

Source as independent or interactive effect

m <- sdmTMB(temp ~ area*year_f + source_f + s(yday, by = area, bs = "cc"), 
            data = df,
            family = student(df = 5),
            spatial = "off",
            spatiotemporal = "off",
            knots = list(yday = c(0.5, 364.5)),
            control = sdmTMBcontrol(newton_loops = 1))

Try alternative degrees of freedom (not evaluated)

#df=20 (effectively gaussian)
m2 <- sdmTMB(temp ~ area*year_f + source_f + s(yday, by = area, bs = "cc"), 
            data = d,
            family = student(df = 20),
            spatial = "off",
            spatiotemporal = "off",
            knots = list(yday = c(0.5, 364.5)),
            control = sdmTMBcontrol(newton_loops = 1))

# df = 9
m3 <- sdmTMB(temp ~ area*year_f + source_f + s(yday, by = area, bs = "cc"), 
            data = d,
            family = student(df = 9),
            spatial = "off",
            spatiotemporal = "off",
            knots = list(yday = c(0.5, 364.5)),
            control = sdmTMBcontrol(newton_loops = 1))

# df=4
m4 <- sdmTMB(temp ~ area*year_f + source_f + s(yday, by = area, bs = "cc"), 
            data = d,
            family = student(df = 4),
            spatial = "off",
            spatiotemporal = "off",
            knots = list(yday = c(0.5, 364.5)),
            control = sdmTMBcontrol(newton_loops = 1))

# Plot all residuals
mcmc_res <- residuals(m, type = "mle-mcmc",
                      mcmc_samples = sdmTMBextra::predict_mle_mcmc(m,
                                                                   mcmc_iter = 201,
                                                                   mcmc_warmup = 200))

mcmc_res20 <- residuals(m2, type = "mle-mcmc",
                        mcmc_samples = sdmTMBextra::predict_mle_mcmc(m2,
                                                                     mcmc_iter = 201,
                                                                     mcmc_warmup = 200))

mcmc_res9 <- residuals(m3, type = "mle-mcmc",
                       mcmc_samples = sdmTMBextra::predict_mle_mcmc(m3,
                                                                    mcmc_iter = 201,
                                                                    mcmc_warmup = 200))

mcmc_res4 <- residuals(m4, type = "mle-mcmc",
                       mcmc_samples = sdmTMBextra::predict_mle_mcmc(m4,
                                                                    mcmc_iter = 201,
                                                                    mcmc_warmup = 200))

dres <- df |> mutate("df=5" = mcmc_res,
                     "df=20" = mcmc_res20,
                     "df=9" = mcmc_res9,
                     "df=4" = mcmc_res4) |> 
  select(`df=5`, `df=20`, `df=9`, `df=4`) |> 
  pivot_longer(everything())

ggplot(dres, aes(sample = value)) +
  stat_qq(size = 0.75, shape = 21, fill = NA) +
  facet_wrap(~factor(name, levels = c("df=20", "df=9", "df=5", "df=4"))) +
  stat_qq_line() +
  labs(y = "Sample Quantiles", x = "Theoretical Quantiles") + 
  theme(aspect.ratio = 1)

Check fit

sanity(m)
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
mcmc_res <- residuals(m, type = "mle-mcmc",
                      mcmc_samples = sdmTMBextra::predict_mle_mcmc(m,
                                                                   mcmc_iter = 201,
                                                                   mcmc_warmup = 200))

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.013802 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 138.02 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 157.592 seconds (Warm-up)
Chain 1:                0.151561 seconds (Sampling)
Chain 1:                157.744 seconds (Total)
Chain 1: 
ggplot(df, aes(sample = mcmc_res)) +
  stat_qq() +
  stat_qq_line() +
  labs(y = "Sample Quantiles", x = "Theoretical Quantiles") + 
  theme(aspect.ratio = 1)

ggsave(paste0(home, "/figures/supp/full_model/qq_temp.pdf"), width = 11, height = 11, units = "cm")

summary(m)
Model fit by ML ['sdmTMB']
Formula: temp ~ area * year_f + source_f + s(yday, by = area, bs = "cc")
Mesh: NULL (isotropic covariance)
Data: df
Family: student(link = 'identity')
 
                     coef.est coef.se
(Intercept)              8.83    0.46
areaFB                   0.06    0.73
areaFM                   0.12    0.77
areaHO                  -0.22    0.80
areaJM                   0.67    0.75
areaMU                   1.29    0.72
areaRA                  -0.78    0.79
areaSI_EK                1.94    0.73
areaTH                   2.01    0.67
areaVN                   1.43    0.66
year_f1941              -0.28    0.69
year_f1942              -0.29    0.66
year_f1943               1.34    0.67
year_f1944               1.35    0.67
year_f1945               1.69    0.66
year_f1946               1.92    0.69
year_f1947               0.86    0.68
year_f1948               1.40    0.66
year_f1949               1.38    0.66
year_f1950               1.16    0.65
year_f1951               1.05    0.67
year_f1952               0.22    0.66
year_f1953               1.62    0.66
year_f1954               1.62    0.66
year_f1955               1.57    0.67
year_f1956               1.07    0.67
year_f1957               1.25    0.67
year_f1958               1.14    0.68
year_f1959               2.26    0.67
year_f1960               1.27    0.65
year_f1961               2.36    0.68
year_f1962               1.42    0.69
year_f1963               1.58    0.66
year_f1964               1.33    0.67
year_f1965               1.26    0.67
year_f1966               1.08    0.65
year_f1967               1.93    0.66
year_f1968               1.63    0.65
year_f1969               1.27    0.65
year_f1970               1.25    0.65
year_f1971               1.19    0.65
year_f1972               1.43    0.66
year_f1973               2.03    0.70
year_f1974               1.86    0.65
year_f1975               2.23    0.65
year_f1976               0.60    0.65
year_f1977               0.62    0.65
year_f1978               0.47    0.65
year_f1979               0.10    0.66
year_f1980               0.57    0.66
year_f1981               0.54    0.65
year_f1982               0.73    0.66
year_f1983               1.30    0.65
year_f1984               1.29    0.65
year_f1985              -0.11    0.65
year_f1986               0.39    0.66
year_f1987              -0.82    0.66
year_f1988               1.34    0.67
year_f1989               2.10    0.66
year_f1990               2.33    0.68
year_f1991               1.56    0.53
year_f1992               1.10    0.51
year_f1993               0.86    0.50
year_f1994               1.33    0.50
year_f1995               2.27    0.51
year_f1996               0.71    0.51
year_f1997               4.65    0.51
year_f1998               0.71    0.48
year_f1999               2.31    0.48
year_f2000               2.06    0.48
year_f2001               2.82    0.48
year_f2002               4.05    0.48
year_f2003               3.30    0.49
year_f2004               2.18    0.48
year_f2005               2.85    0.48
year_f2006               0.95    0.52
year_f2007               1.36    0.52
year_f2008               2.22    0.52
year_f2009               2.53    0.52
year_f2010               1.33    0.67
year_f2011               1.75    0.68
year_f2012               1.74    0.65
year_f2013               2.13    0.66
year_f2014               2.50    0.65
year_f2015               2.55    0.65
year_f2016               2.65    0.65
year_f2017               2.20    0.65
year_f2018               2.80    0.66
year_f2019               2.14    0.65
year_f2020               3.22    0.64
year_f2021               2.34    0.66
year_f2022               2.28    1.20
source_ffishing          0.32    0.04
source_flogger           0.15    0.04
areaFB:year_f1941        0.06    1.11
areaFM:year_f1941        0.07    1.17
areaHO:year_f1941        0.41    1.18
areaJM:year_f1941       -0.18    1.15
areaMU:year_f1941       -0.16    1.07
areaRA:year_f1941        0.09    1.14
areaSI_EK:year_f1941    -0.44    1.09
areaTH:year_f1941       -0.59    0.98
areaVN:year_f1941       -0.52    0.96
areaFB:year_f1942        0.14    1.04
areaFM:year_f1942        0.16    1.10
areaHO:year_f1942        0.38    1.12
areaJM:year_f1942       -0.04    1.09
areaMU:year_f1942       -0.19    1.01
areaRA:year_f1942        0.05    1.09
areaSI_EK:year_f1942    -0.18    1.07
areaTH:year_f1942        0.08    0.98
areaVN:year_f1942       -0.24    0.96
areaFB:year_f1943       -0.58    1.05
areaFM:year_f1943       -0.79    1.09
areaHO:year_f1943       -0.87    1.13
areaJM:year_f1943       -0.64    1.07
areaMU:year_f1943       -0.40    1.01
areaRA:year_f1943       -0.87    1.13
areaSI_EK:year_f1943    -0.31    1.03
areaTH:year_f1943        0.05    0.95
areaVN:year_f1943       -0.04    0.95
areaFB:year_f1944       -0.32    1.08
areaFM:year_f1944       -0.31    1.15
areaHO:year_f1944       -0.44    1.18
areaJM:year_f1944       -0.41    1.11
areaMU:year_f1944       -0.24    1.02
areaRA:year_f1944       -0.80    1.14
areaSI_EK:year_f1944    -0.15    1.04
areaTH:year_f1944       -0.20    0.95
areaVN:year_f1944       -0.21    0.94
areaFB:year_f1945       -0.23    1.03
areaFM:year_f1945       -0.43    1.08
areaHO:year_f1945       -0.75    1.12
areaJM:year_f1945       -0.18    1.06
areaMU:year_f1945       -0.16    0.99
areaRA:year_f1945       -1.16    1.11
areaSI_EK:year_f1945     0.07    1.02
areaTH:year_f1945        0.27    0.94
areaVN:year_f1945        0.21    0.93
areaFB:year_f1946        0.13    1.04
areaFM:year_f1946       -0.07    1.10
areaHO:year_f1946       -0.59    1.13
areaJM:year_f1946        0.19    1.07
areaMU:year_f1946       -0.08    1.00
areaRA:year_f1946       -1.68    1.12
areaSI_EK:year_f1946    -0.01    1.03
areaTH:year_f1946       -0.06    0.96
areaVN:year_f1946        0.16    0.96
areaFB:year_f1947       -0.10    1.04
areaFM:year_f1947       -0.24    1.09
areaHO:year_f1947       -1.01    1.13
areaJM:year_f1947        0.03    1.08
areaMU:year_f1947       -0.02    1.01
areaRA:year_f1947       -1.87    1.18
areaSI_EK:year_f1947     0.21    1.06
areaTH:year_f1947        0.22    0.99
areaVN:year_f1947        0.37    0.98
areaFB:year_f1948       -0.11    1.03
areaFM:year_f1948       -0.15    1.09
areaHO:year_f1948       -0.23    1.12
areaJM:year_f1948       -0.13    1.06
areaMU:year_f1948       -0.27    1.00
areaRA:year_f1948       -0.85    1.09
areaSI_EK:year_f1948    -0.26    1.03
areaTH:year_f1948       -0.28    0.94
areaVN:year_f1948       -0.04    0.93
areaFB:year_f1949       -0.20    1.08
areaFM:year_f1949       -0.40    1.14
areaHO:year_f1949       -0.54    1.20
areaJM:year_f1949       -0.18    1.11
areaMU:year_f1949        0.06    1.03
areaRA:year_f1949       -1.37    1.17
areaSI_EK:year_f1949     0.32    1.08
areaTH:year_f1949        0.44    0.98
areaVN:year_f1949        0.39    0.94
areaFB:year_f1950        0.01    1.04
areaFM:year_f1950       -0.12    1.11
areaHO:year_f1950       -0.22    1.13
areaJM:year_f1950        0.22    1.07
areaMU:year_f1950        0.01    0.98
areaRA:year_f1950       -0.61    1.14
areaSI_EK:year_f1950     0.24    1.01
areaTH:year_f1950        0.28    0.93
areaVN:year_f1950        0.37    0.92
areaFB:year_f1951       -0.28    1.11
areaFM:year_f1951       -0.40    1.19
areaHO:year_f1951       -0.43    1.21
areaJM:year_f1951       -0.30    1.14
areaMU:year_f1951       -0.03    1.05
areaRA:year_f1951       -1.12    1.17
areaSI_EK:year_f1951     0.12    1.09
areaTH:year_f1951        0.00    0.98
areaVN:year_f1951        0.10    0.95
areaFB:year_f1952       -0.11    1.02
areaFM:year_f1952       -0.23    1.06
areaHO:year_f1952       -0.40    1.13
areaJM:year_f1952        0.15    1.05
areaMU:year_f1952       -0.10    1.00
areaRA:year_f1952       -0.95    1.11
areaSI_EK:year_f1952     0.07    1.03
areaTH:year_f1952        0.31    0.94
areaVN:year_f1952        0.37    0.94
areaFB:year_f1953        0.15    1.04
areaFM:year_f1953        0.06    1.07
areaHO:year_f1953       -0.31    1.14
areaJM:year_f1953       -0.02    1.07
areaMU:year_f1953       -0.12    1.00
areaRA:year_f1953       -1.17    1.11
areaSI_EK:year_f1953    -0.31    1.02
areaTH:year_f1953       -0.16    0.95
areaVN:year_f1953       -0.07    0.94
areaFB:year_f1954        0.17    1.06
areaFM:year_f1954        0.08    1.10
areaHO:year_f1954       -0.16    1.15
areaJM:year_f1954       -0.10    1.08
areaMU:year_f1954       -0.40    1.03
areaRA:year_f1954       -1.37    1.09
areaSI_EK:year_f1954    -0.68    1.05
areaTH:year_f1954       -0.60    0.96
areaVN:year_f1954       -0.43    0.94
areaFB:year_f1955        0.71    1.07
areaFM:year_f1955        0.81    1.10
areaHO:year_f1955        0.26    1.16
areaJM:year_f1955        0.65    1.10
areaMU:year_f1955        0.11    1.04
areaRA:year_f1955       -1.52    1.10
areaSI_EK:year_f1955     0.08    1.07
areaTH:year_f1955       -0.25    0.99
areaVN:year_f1955       -0.35    0.96
areaFB:year_f1956        0.04    1.07
areaFM:year_f1956       -0.10    1.11
areaHO:year_f1956       -0.05    1.19
areaJM:year_f1956       -0.30    1.09
areaMU:year_f1956       -0.64    1.04
areaRA:year_f1956       -1.07    1.11
areaSI_EK:year_f1956    -0.94    1.05
areaTH:year_f1956       -0.78    0.97
areaVN:year_f1956       -0.70    0.94
areaFB:year_f1957        0.41    1.07
areaFM:year_f1957        0.42    1.10
areaHO:year_f1957        0.34    1.14
areaJM:year_f1957        0.30    1.09
areaMU:year_f1957       -0.26    1.05
areaRA:year_f1957       -0.92    1.11
areaSI_EK:year_f1957    -0.35    1.06
areaTH:year_f1957       -0.12    0.96
areaVN:year_f1957       -0.12    0.94
areaFB:year_f1958        0.46    1.14
areaFM:year_f1958        0.42    1.18
areaHO:year_f1958        0.51    1.19
areaJM:year_f1958        0.18    1.18
areaMU:year_f1958       -0.15    1.11
areaRA:year_f1958       -0.82    1.14
areaSI_EK:year_f1958    -0.27    1.15
areaTH:year_f1958       -0.18    1.01
areaVN:year_f1958       -0.26    0.96
areaFB:year_f1959        0.30    1.03
areaFM:year_f1959        0.24    1.07
areaHO:year_f1959       -0.16    1.11
areaJM:year_f1959        0.28    1.06
areaMU:year_f1959       -0.18    1.01
areaRA:year_f1959       -1.66    1.09
areaSI_EK:year_f1959    -0.17    1.03
areaTH:year_f1959       -0.07    0.94
areaVN:year_f1959       -0.06    0.94
areaFB:year_f1960        0.25    1.02
areaFM:year_f1960        0.22    1.06
areaHO:year_f1960       -0.13    1.11
areaJM:year_f1960        0.13    1.05
areaMU:year_f1960       -0.19    1.00
areaRA:year_f1960       -1.04    1.08
areaSI_EK:year_f1960    -0.28    1.02
areaTH:year_f1960       -0.15    0.93
areaVN:year_f1960       -0.12    0.92
areaFB:year_f1961        0.03    1.10
areaFM:year_f1961       -0.16    1.14
areaHO:year_f1961       -0.62    1.17
areaJM:year_f1961       -0.15    1.12
areaMU:year_f1961       -0.24    1.06
areaRA:year_f1961       -1.72    1.13
areaSI_EK:year_f1961    -0.50    1.07
areaTH:year_f1961       -0.43    0.97
areaVN:year_f1961       -0.27    0.95
areaFB:year_f1962        0.18    1.15
areaFM:year_f1962       -0.06    1.19
areaHO:year_f1962       -0.07    1.21
areaJM:year_f1962       -0.16    1.18
areaMU:year_f1962       -0.47    1.12
areaRA:year_f1962       -1.17    1.14
areaSI_EK:year_f1962    -0.96    1.15
areaTH:year_f1962       -0.84    1.01
areaVN:year_f1962       -0.57    0.97
areaFB:year_f1963        0.13    1.03
areaFM:year_f1963        0.01    1.09
areaHO:year_f1963       -0.24    1.14
areaJM:year_f1963       -0.20    1.07
areaMU:year_f1963       -0.30    0.99
areaRA:year_f1963       -1.08    1.12
areaSI_EK:year_f1963    -0.65    1.04
areaTH:year_f1963       -0.89    0.96
areaVN:year_f1963       -0.58    0.94
areaFB:year_f1964        0.06    1.09
areaFM:year_f1964       -0.18    1.13
areaHO:year_f1964       -0.05    1.17
areaJM:year_f1964       -0.28    1.10
areaMU:year_f1964       -0.42    1.05
areaRA:year_f1964       -0.93    1.14
areaSI_EK:year_f1964    -0.73    1.07
areaTH:year_f1964       -0.35    0.97
areaVN:year_f1964       -0.39    0.94
areaFB:year_f1965        0.23    1.09
areaFM:year_f1965        0.07    1.14
areaHO:year_f1965        0.22    1.17
areaJM:year_f1965        0.04    1.13
areaMU:year_f1965       -0.50    1.06
areaRA:year_f1965       -0.88    1.12
areaSI_EK:year_f1965    -0.86    1.09
areaTH:year_f1965       -0.70    0.98
areaVN:year_f1965       -0.62    0.95
areaFB:year_f1966        0.09    1.03
areaFM:year_f1966        0.04    1.06
areaHO:year_f1966       -0.37    1.13
areaJM:year_f1966       -0.05    1.05
areaMU:year_f1966       -0.11    1.00
areaRA:year_f1966       -1.25    1.08
areaSI_EK:year_f1966    -0.19    1.02
areaTH:year_f1966       -0.10    0.95
areaVN:year_f1966       -0.10    0.94
areaFB:year_f1967        0.26    1.07
areaFM:year_f1967        0.16    1.10
areaHO:year_f1967        0.03    1.15
areaJM:year_f1967        0.08    1.11
areaMU:year_f1967       -0.11    1.02
areaRA:year_f1967       -1.08    1.13
areaSI_EK:year_f1967    -0.23    1.05
areaTH:year_f1967       -0.16    0.95
areaVN:year_f1967       -0.15    0.93
areaFB:year_f1968        0.28    1.04
areaFM:year_f1968        0.20    1.08
areaHO:year_f1968        0.15    1.14
areaJM:year_f1968        0.15    1.08
areaMU:year_f1968       -0.21    1.00
areaRA:year_f1968       -1.21    1.12
areaSI_EK:year_f1968    -0.33    1.03
areaTH:year_f1968       -0.19    0.94
areaVN:year_f1968       -0.21    0.92
areaFB:year_f1969        0.40    1.02
areaFM:year_f1969        0.40    1.07
areaHO:year_f1969        0.02    1.11
areaJM:year_f1969        0.44    1.06
areaMU:year_f1969        0.06    0.99
areaRA:year_f1969       -1.04    1.09
areaSI_EK:year_f1969     0.24    1.03
areaTH:year_f1969        0.24    0.94
areaVN:year_f1969        0.14    0.93
areaFB:year_f1970        0.11    1.02
areaFM:year_f1970        0.03    1.06
areaHO:year_f1970       -0.10    1.12
areaJM:year_f1970       -0.02    1.05
areaMU:year_f1970       -0.22    0.99
areaRA:year_f1970       -0.86    1.09
areaSI_EK:year_f1970    -0.28    1.01
areaTH:year_f1970       -0.27    0.94
areaVN:year_f1970       -0.20    0.92
areaFB:year_f1971        0.23    1.05
areaFM:year_f1971        0.20    1.09
areaHO:year_f1971       -0.25    1.14
areaJM:year_f1971        0.30    1.08
areaMU:year_f1971       -0.04    1.02
areaRA:year_f1971       -1.26    1.08
areaSI_EK:year_f1971     0.05    1.04
areaTH:year_f1971        0.19    0.94
areaVN:year_f1971        0.19    0.92
areaFB:year_f1972        0.50    1.04
areaFM:year_f1972        0.53    1.10
areaHO:year_f1972       -0.32    1.12
areaJM:year_f1972        0.39    0.94
areaMU:year_f1972        0.39    1.00
areaRA:year_f1972       -1.02    1.10
areaSI_EK:year_f1972     0.41    1.03
areaTH:year_f1972        0.26    0.94
areaVN:year_f1972        0.38    0.94
areaFB:year_f1973        0.22    1.04
areaFM:year_f1973        0.18    1.07
areaHO:year_f1973       -0.01    1.12
areaJM:year_f1973       -2.76    0.96
areaMU:year_f1973       -0.23    1.04
areaRA:year_f1973       -0.95    1.09
areaSI_EK:year_f1973    -0.15    1.04
areaTH:year_f1973       -0.04    0.98
areaVN:year_f1973        0.11    0.99
areaFB:year_f1974        0.34    1.04
areaFM:year_f1974        0.23    1.08
areaHO:year_f1974        0.15    1.14
areaJM:year_f1974       -3.72    0.92
areaMU:year_f1974       -0.06    1.00
areaRA:year_f1974       -1.27    1.09
areaSI_EK:year_f1974    -0.29    1.03
areaTH:year_f1974       -0.39    0.94
areaVN:year_f1974       -0.15    0.92
areaFB:year_f1975        0.51    1.04
areaFM:year_f1975       -0.64    0.98
areaHO:year_f1975        0.52    1.13
areaJM:year_f1975       -1.91    0.93
areaMU:year_f1975        0.02    1.01
areaRA:year_f1975       -0.60    1.09
areaSI_EK:year_f1975     0.11    1.04
areaTH:year_f1975        0.02    0.95
areaVN:year_f1975        0.07    0.93
areaFB:year_f1976       -0.43    0.96
areaFM:year_f1976        0.70    0.97
areaHO:year_f1976        0.54    1.12
areaJM:year_f1976        0.02    0.93
areaMU:year_f1976        0.05    1.01
areaRA:year_f1976       -0.43    1.10
areaSI_EK:year_f1976     0.36    1.02
areaTH:year_f1976        0.63    0.94
areaVN:year_f1976        0.38    0.92
areaFB:year_f1977       -2.00    0.93
areaFM:year_f1977       -0.77    0.95
areaHO:year_f1977       -0.06    1.15
areaJM:year_f1977       -1.52    0.91
areaMU:year_f1977        0.07    1.02
areaRA:year_f1977       -0.71    1.13
areaSI_EK:year_f1977     0.09    1.03
areaTH:year_f1977        0.31    0.94
areaVN:year_f1977        0.23    0.92
areaFB:year_f1978       -1.29    0.91
areaFM:year_f1978       -0.85    0.99
areaHO:year_f1978        0.00    1.13
areaJM:year_f1978       -1.12    0.93
areaMU:year_f1978        0.00    1.02
areaRA:year_f1978       -0.97    1.08
areaSI_EK:year_f1978     0.00    1.04
areaTH:year_f1978        0.27    0.95
areaVN:year_f1978        0.19    0.92
areaFB:year_f1979       -0.18    0.93
areaFM:year_f1979       -0.27    0.98
areaHO:year_f1979       -0.15    1.12
areaJM:year_f1979       -2.74    0.93
areaMU:year_f1979        0.13    0.99
areaRA:year_f1979       -0.57    1.09
areaSI_EK:year_f1979     0.25    1.02
areaTH:year_f1979        0.24    0.95
areaVN:year_f1979        0.32    0.93
areaFB:year_f1980       -1.01    0.97
areaFM:year_f1980        0.41    1.01
areaHO:year_f1980        0.13    1.15
areaJM:year_f1980       -0.25    0.93
areaMU:year_f1980        0.02    0.99
areaRA:year_f1980       -0.63    1.11
areaSI_EK:year_f1980     0.12    1.02
areaTH:year_f1980        0.09    0.94
areaVN:year_f1980        0.15    0.93
areaFB:year_f1981       -0.94    0.93
areaFM:year_f1981        0.23    1.00
areaHO:year_f1981       -0.35    1.14
areaJM:year_f1981       -0.68    0.92
areaMU:year_f1981        0.14    1.00
areaRA:year_f1981       -1.28    1.11
areaSI_EK:year_f1981     0.17    1.03
areaTH:year_f1981        0.23    0.94
areaVN:year_f1981        0.24    0.92
areaFB:year_f1982       -1.49    0.93
areaFM:year_f1982       -0.36    0.98
areaHO:year_f1982       -0.52    1.16
areaJM:year_f1982       -2.42    0.93
areaMU:year_f1982        0.11    1.01
areaRA:year_f1982       -0.67    1.13
areaSI_EK:year_f1982     0.24    1.04
areaTH:year_f1982        0.49    0.96
areaVN:year_f1982        0.43    0.94
areaFB:year_f1983       -1.37    0.91
areaFM:year_f1983        0.11    0.92
areaHO:year_f1983       -0.14    1.11
areaJM:year_f1983       -2.30    0.92
areaMU:year_f1983        0.01    1.00
areaRA:year_f1983       -0.96    1.07
areaSI_EK:year_f1983     0.10    1.02
areaTH:year_f1983        0.16    0.94
areaVN:year_f1983        0.18    0.93
areaFB:year_f1984       -0.05    0.95
areaFM:year_f1984        1.31    0.92
areaHO:year_f1984       -0.36    1.11
areaJM:year_f1984       -0.63    0.91
areaMU:year_f1984        0.02    0.99
areaRA:year_f1984       -1.28    1.08
areaSI_EK:year_f1984     0.03    1.04
areaTH:year_f1984       -0.07    0.95
areaVN:year_f1984        0.08    0.92
areaFB:year_f1985        0.59    0.90
areaFM:year_f1985        3.92    0.92
areaHO:year_f1985        0.06    1.12
areaJM:year_f1985       -4.40    0.92
areaMU:year_f1985        0.11    1.00
areaRA:year_f1985       -0.44    1.11
areaSI_EK:year_f1985     0.29    1.02
areaTH:year_f1985        0.21    0.94
areaVN:year_f1985        0.30    0.92
areaFB:year_f1986       -1.36    0.92
areaFM:year_f1986        9.74    0.92
areaHO:year_f1986       -0.21    1.09
areaJM:year_f1986       -3.27    0.97
areaMU:year_f1986       -0.08    1.00
areaRA:year_f1986       -0.24    1.07
areaSI_EK:year_f1986    -0.17    1.03
areaTH:year_f1986       -0.16    0.95
areaVN:year_f1986        0.01    0.93
areaFB:year_f1987       -1.05    0.92
areaFM:year_f1987        5.26    0.93
areaHO:year_f1987       -0.20    1.20
areaJM:year_f1987       -0.38    0.91
areaMU:year_f1987        0.06    1.03
areaRA:year_f1987       -0.05    1.17
areaSI_EK:year_f1987    -1.53    0.92
areaTH:year_f1987        0.54    0.97
areaVN:year_f1987        0.51    0.94
areaFB:year_f1988       -2.14    0.90
areaFM:year_f1988        6.88    0.93
areaHO:year_f1988        0.18    1.11
areaJM:year_f1988       -2.41    0.91
areaMU:year_f1988       -0.08    1.01
areaRA:year_f1988       -1.23    1.09
areaSI_EK:year_f1988    -3.87    0.94
areaTH:year_f1988       -0.09    0.95
areaVN:year_f1988       -0.09    0.94
areaFB:year_f1989       -3.13    0.90
areaFM:year_f1989        0.01    0.91
areaHO:year_f1989       -1.62    0.95
areaJM:year_f1989       -1.47    0.89
areaMU:year_f1989       -0.19    0.99
areaRA:year_f1989       -1.18    1.08
areaSI_EK:year_f1989    -1.38    0.88
areaTH:year_f1989       -0.01    0.94
areaVN:year_f1989        0.00    0.93
areaFB:year_f1990       -2.11    0.90
areaFM:year_f1990        2.69    0.92
areaHO:year_f1990       -1.20    0.95
areaJM:year_f1990        1.27    0.91
areaMU:year_f1990       -0.26    1.01
areaRA:year_f1990       -1.06    1.08
areaSI_EK:year_f1990    -0.92    0.89
areaTH:year_f1990       -0.07    0.95
areaVN:year_f1990        0.00    0.96
areaFB:year_f1991        0.60    0.79
areaFM:year_f1991        3.52    0.84
areaHO:year_f1991       -0.35    0.85
areaJM:year_f1991       -0.27    0.81
areaMU:year_f1991       -0.01    0.80
areaRA:year_f1991        0.39    0.85
areaSI_EK:year_f1991    -0.61    0.85
areaTH:year_f1991       -0.15    0.87
areaVN:year_f1991       -0.01    0.85
areaFB:year_f1992        0.90    0.77
areaFM:year_f1992        1.29    0.80
areaHO:year_f1992        0.04    0.84
areaJM:year_f1992        1.37    0.84
areaMU:year_f1992        0.99    0.82
areaRA:year_f1992        1.59    0.87
areaSI_EK:year_f1992    -1.07    0.77
areaTH:year_f1992        0.63    0.85
areaVN:year_f1992        0.75    0.85
areaFB:year_f1993        0.59    0.76
areaFM:year_f1993        0.73    0.80
areaHO:year_f1993       -0.80    0.83
areaJM:year_f1993        0.81    0.78
areaMU:year_f1993       -0.12    0.82
areaRA:year_f1993       -0.53    0.99
areaSI_EK:year_f1993    -0.44    0.77
areaTH:year_f1993        0.02    0.84
areaVN:year_f1993        0.10    0.84
areaFB:year_f1994        0.89    0.76
areaFM:year_f1994        2.49    0.79
areaHO:year_f1994        1.06    0.83
areaJM:year_f1994        1.39    0.78
areaMU:year_f1994       -1.41    0.83
areaRA:year_f1994        0.07    0.84
areaSI_EK:year_f1994     0.60    0.77
areaTH:year_f1994        0.08    0.84
areaVN:year_f1994        0.04    0.83
areaFB:year_f1995       -0.18    0.77
areaFM:year_f1995        0.86    0.81
areaHO:year_f1995       -1.56    0.84
areaJM:year_f1995        0.93    0.79
areaMU:year_f1995       -0.34    0.82
areaRA:year_f1995       -4.62    0.86
areaSI_EK:year_f1995    -3.40    0.79
areaTH:year_f1995       -0.51    0.85
areaVN:year_f1995       -0.11    0.75
areaFB:year_f1996        1.27    0.77
areaFM:year_f1996        2.84    0.80
areaHO:year_f1996        0.22    0.84
areaJM:year_f1996        0.69    0.82
areaMU:year_f1996       -0.90    0.76
areaRA:year_f1996        0.00    0.85
areaSI_EK:year_f1996    -1.40    0.77
areaTH:year_f1996       -0.13    0.87
areaVN:year_f1996       -2.28    0.77
areaFB:year_f1997       -1.57    0.78
areaFM:year_f1997       -1.83    0.81
areaHO:year_f1997       -2.85    0.90
areaJM:year_f1997       -0.74    0.83
areaMU:year_f1997       -3.50    0.76
areaRA:year_f1997       -2.07    0.87
areaSI_EK:year_f1997    -2.78    0.77
areaTH:year_f1997       -2.83    0.86
areaVN:year_f1997        0.19    0.73
areaFB:year_f1998        1.02    0.75
areaFM:year_f1998        1.26    0.81
areaHO:year_f1998       -0.21    0.82
areaJM:year_f1998        1.06    0.77
areaMU:year_f1998       -0.91    0.73
areaRA:year_f1998       -0.05    0.81
areaSI_EK:year_f1998    -1.20    0.76
areaTH:year_f1998        0.51    0.83
areaVN:year_f1998       -6.00    0.72
areaFB:year_f1999        1.28    0.75
areaFM:year_f1999        0.65    0.80
areaHO:year_f1999       -3.05    0.88
areaJM:year_f1999        1.10    0.77
areaMU:year_f1999       -1.38    0.73
areaRA:year_f1999       -1.62    0.81
areaSI_EK:year_f1999    -2.08    0.75
areaTH:year_f1999       -0.10    0.83
areaVN:year_f1999       -2.14    0.72
areaFB:year_f2000        1.23    0.75
areaFM:year_f2000        0.99    0.78
areaHO:year_f2000       -0.29    0.82
areaJM:year_f2000        1.72    0.77
areaMU:year_f2000       -0.39    0.73
areaRA:year_f2000       -0.76    0.80
areaSI_EK:year_f2000    -0.62    0.75
areaTH:year_f2000       -0.04    0.84
areaVN:year_f2000        0.00    0.70
areaFB:year_f2001        0.71    0.75
areaFM:year_f2001       -0.11    0.78
areaHO:year_f2001       -1.34    0.82
areaJM:year_f2001        0.43    0.76
areaMU:year_f2001       -2.00    0.73
areaRA:year_f2001       -0.23    0.80
areaSI_EK:year_f2001    -2.79    0.75
areaTH:year_f2001       -0.85    0.83
areaVN:year_f2001       -4.90    0.72
areaFB:year_f2002        0.61    0.76
areaFM:year_f2002        1.55    0.79
areaHO:year_f2002       -1.78    0.82
areaJM:year_f2002        0.64    0.77
areaMU:year_f2002       -2.88    0.73
areaRA:year_f2002       -1.04    0.82
areaSI_EK:year_f2002    -1.85    0.76
areaTH:year_f2002       -0.12    0.80
areaVN:year_f2002        0.04    0.71
areaFB:year_f2003        0.06    0.76
areaFM:year_f2003        0.72    0.80
areaHO:year_f2003       -0.73    0.83
areaJM:year_f2003       -0.19    0.77
areaMU:year_f2003       -2.71    0.74
areaRA:year_f2003       -1.33    0.82
areaSI_EK:year_f2003    -2.46    0.76
areaTH:year_f2003       -1.52    0.78
areaVN:year_f2003        0.98    0.71
areaFB:year_f2004        0.91    0.75
areaFM:year_f2004        0.83    0.79
areaHO:year_f2004       -0.31    0.82
areaJM:year_f2004        0.95    0.77
areaMU:year_f2004       -1.48    0.73
areaRA:year_f2004       -1.13    0.81
areaSI_EK:year_f2004    -1.26    0.75
areaTH:year_f2004        0.66    0.78
areaVN:year_f2004        1.95    0.72
areaFB:year_f2005        0.63    0.75
areaFM:year_f2005        0.09    0.79
areaHO:year_f2005       -1.09    0.82
areaJM:year_f2005        0.34    0.77
areaMU:year_f2005       -2.06    0.73
areaRA:year_f2005       -0.79    0.81
areaSI_EK:year_f2005    -2.32    0.78
areaTH:year_f2005       -0.68    0.77
areaVN:year_f2005       -1.77    0.70
areaFB:year_f2006        3.37    0.79
areaFM:year_f2006        1.47    0.82
areaHO:year_f2006        1.82    0.85
areaJM:year_f2006        3.91    0.80
areaMU:year_f2006        0.49    0.76
areaRA:year_f2006        3.30    0.95
areaSI_EK:year_f2006     1.68    0.78
areaTH:year_f2006        1.04    0.72
areaVN:year_f2006        2.28    0.74
areaFB:year_f2007        1.73    0.78
areaFM:year_f2007        1.09    0.81
areaHO:year_f2007        0.38    0.84
areaJM:year_f2007        1.85    0.79
areaMU:year_f2007       -0.10    0.76
areaRA:year_f2007        0.07    0.83
areaSI_EK:year_f2007    -0.06    0.77
areaTH:year_f2007        1.23    0.79
areaVN:year_f2007        0.91    0.77
areaFB:year_f2008        1.16    0.78
areaFM:year_f2008        0.88    0.81
areaHO:year_f2008       -0.76    0.84
areaJM:year_f2008        1.17    0.79
areaMU:year_f2008       -0.35    0.76
areaRA:year_f2008       -1.05    0.84
areaSI_EK:year_f2008    -0.26    0.77
areaTH:year_f2008       -0.57    0.80
areaVN:year_f2008       -0.72    0.76
areaFB:year_f2009        0.81    0.78
areaFM:year_f2009       -0.04    0.81
areaHO:year_f2009       -0.77    0.85
areaJM:year_f2009        0.12    0.79
areaMU:year_f2009       -1.41    0.76
areaRA:year_f2009        0.16    0.84
areaSI_EK:year_f2009    -2.41    0.78
areaTH:year_f2009       -0.57    0.79
areaVN:year_f2009        0.46    0.79
areaFB:year_f2010        1.30    0.89
areaFM:year_f2010        0.83    0.92
areaHO:year_f2010       -0.06    0.95
areaJM:year_f2010        1.55    0.90
areaMU:year_f2010       -0.55    0.87
areaRA:year_f2010        0.87    0.94
areaSI_EK:year_f2010     0.16    0.96
areaTH:year_f2010       -2.19    0.86
areaVN:year_f2010        2.00    0.92
areaFB:year_f2011        2.13    0.89
areaFM:year_f2011        1.67    0.92
areaHO:year_f2011        1.07    0.95
areaJM:year_f2011        1.67    0.90
areaMU:year_f2011       -0.23    0.88
areaRA:year_f2011        1.36    0.94
areaSI_EK:year_f2011    -0.58    0.89
areaTH:year_f2011       -1.62    0.84
areaVN:year_f2011        0.60    0.89
areaFB:year_f2012        1.25    0.87
areaFM:year_f2012        1.11    0.90
areaHO:year_f2012       -0.15    0.93
areaJM:year_f2012        0.85    0.88
areaMU:year_f2012       -0.69    0.85
areaRA:year_f2012        0.18    0.97
areaSI_EK:year_f2012    -1.06    0.87
areaTH:year_f2012       -0.87    0.81
areaVN:year_f2012       -2.31    0.84
areaFB:year_f2013        1.78    0.88
areaFM:year_f2013        1.59    0.91
areaHO:year_f2013        0.58    0.94
areaJM:year_f2013        1.72    0.89
areaMU:year_f2013       -0.10    0.87
areaRA:year_f2013       -0.15    0.98
areaSI_EK:year_f2013    -0.90    0.88
areaTH:year_f2013       -0.24    0.83
areaVN:year_f2013        1.44    0.84
areaFB:year_f2014        1.66    0.88
areaFM:year_f2014        0.98    0.90
areaHO:year_f2014        0.15    0.93
areaJM:year_f2014        2.15    0.89
areaMU:year_f2014        0.23    0.86
areaRA:year_f2014        1.10    0.93
areaSI_EK:year_f2014    -0.15    0.87
areaTH:year_f2014        0.21    0.82
areaVN:year_f2014        1.57    0.90
areaFB:year_f2015        0.97    0.87
areaFM:year_f2015        0.50    0.90
areaHO:year_f2015       -0.55    0.94
areaJM:year_f2015        1.30    0.88
areaMU:year_f2015       -0.19    0.86
areaRA:year_f2015       -0.51    0.93
areaSI_EK:year_f2015    -0.62    0.87
areaTH:year_f2015       -0.13    0.82
areaVN:year_f2015        0.00    0.88
areaFB:year_f2016        1.03    0.88
areaFM:year_f2016       -0.29    0.91
areaHO:year_f2016       -0.85    0.94
areaJM:year_f2016        0.85    0.89
areaMU:year_f2016       -0.93    0.86
areaRA:year_f2016        0.15    0.93
areaSI_EK:year_f2016    -1.11    0.87
areaTH:year_f2016       -0.52    0.82
areaVN:year_f2016       -2.56    1.12
areaFB:year_f2017        0.74    0.87
areaFM:year_f2017        0.24    0.90
areaHO:year_f2017       -0.80    0.93
areaJM:year_f2017        1.08    0.88
areaMU:year_f2017       -0.52    0.85
areaRA:year_f2017       -0.75    0.92
areaSI_EK:year_f2017    -0.78    0.86
areaTH:year_f2017       -1.06    0.81
areaVN:year_f2017       -0.87    0.88
areaFB:year_f2018        1.43    0.88
areaFM:year_f2018        1.55    0.91
areaHO:year_f2018       -0.84    0.93
areaJM:year_f2018        1.33    0.89
areaMU:year_f2018       -1.38    0.86
areaRA:year_f2018       -0.20    0.93
areaSI_EK:year_f2018    -1.55    0.87
areaTH:year_f2018        0.09    0.88
areaVN:year_f2018       -0.59    0.89
areaFB:year_f2019       -0.04    0.95
areaFM:year_f2019        0.85    0.90
areaHO:year_f2019       -0.55    0.94
areaJM:year_f2019        1.31    0.89
areaMU:year_f2019       -0.75    0.86
areaRA:year_f2019       -0.30    0.93
areaSI_EK:year_f2019    -0.28    0.87
areaTH:year_f2019       -0.57    0.82
areaVN:year_f2019        0.05    0.92
areaFB:year_f2020       -0.11    0.93
areaFM:year_f2020       -0.60    0.91
areaHO:year_f2020       -0.74    1.07
areaJM:year_f2020        1.05    0.92
areaMU:year_f2020        0.46    0.87
areaRA:year_f2020       -0.38    0.98
areaSI_EK:year_f2020    -1.67    1.04
areaTH:year_f2020        1.39    0.92
areaVN:year_f2020       -0.20    0.86
areaFB:year_f2021        0.24    1.02
areaFM:year_f2021        0.66    0.97
areaHO:year_f2021       -0.92    1.04
areaJM:year_f2021        1.33    0.92
areaMU:year_f2021        1.63    0.89
areaRA:year_f2021       -2.03    0.99
areaSI_EK:year_f2021    -0.54    1.03
areaTH:year_f2021       -0.91    0.91
areaVN:year_f2021       -0.04    0.92
areaFB:year_f2022        0.49    1.72
areaFM:year_f2022        0.61    1.40
areaHO:year_f2022       -0.91    1.48
areaJM:year_f2022        0.44    1.41
areaMU:year_f2022        1.78    1.33
areaRA:year_f2022       -1.18    1.43
areaSI_EK:year_f2022     0.56    1.54
areaTH:year_f2022        1.38    1.37
areaVN:year_f2022       -4.31    1.36

Smooth terms:
                     Std. Dev.
sds(yday):areaBS)         0.90
sds(yday):areaFB)         1.04
sds(yday):areaFM)         0.95
sds(yday):areaHO)         1.17
sds(yday):areaJM)         1.30
sds(yday):areaMU)         0.87
sds(yday):areaRA)         1.43
sds(yday):areaSI_EK)      0.75
sds(yday):areaTH)         0.79
sds(yday):areaVN)         0.79

Dispersion parameter: 1.67
ML criterion at convergence: 158107.311

See ?tidy.sdmTMB to extract these values as a data frame.

Predict

# Make a new data frame and predict!
nd <- data.frame(expand.grid(yday = seq(min(df$yday), max(df$yday), by = 1),
                             area = unique(df$area),
                             source = unique(df$source_f),
                             year = unique(df$year))) |>
  mutate(source_f = as.factor(source),
         year_f = as.factor(year)) |> # Left join in growth data column
  left_join(gdat, by = "area") |> 
  mutate(area = as.factor(area),
         growth_dat = ifelse(year >= min & year <= max, "Y", "N"))
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x         0
           > rows only in y  (      2)
           > matched rows     911,340
           >                 =========
           > rows total       911,340
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
# Predict
nd$pred <- predict(m, newdata = nd)$est

In order to have a lower temperature in before-nuclear times (without any data to inform that), we can use the nearby areas.. so FM informs BT prior to nuclear

nd_sub <- nd |> 
  mutate(keep = "N",
         keep = ifelse(area == "FM" & year <= 1980, "Y", keep), # use FM instead of BT
         keep = ifelse(area == "SI_EK" & year <= 1972, "Y", keep)) |> # use SI_EK instead of SI_HA
  filter(keep == "Y") |> # Now change the labels to BT and SI_EK...
  mutate(area = ifelse(area == "FM", "BT", "SI_HA"))
mutate: new variable 'keep' (character) with 2 unique values and 0% NA
filter: removed 830,088 rows (91%), 81,252 rows remaining
mutate: converted 'area' from factor to character (0 new NA)
# Bind rows and plot only the temperature series we will use for growth modelling
nd <- bind_rows(nd, nd_sub) |>
  select(-keep) |> 
  mutate(growth_dat = ifelse(area == "SI_HA" & year %in% c(1966, 1967), "Y", growth_dat)) # SI_EK and SI_HA do not have the same starting years, so we can't use allo columns from SI_EK
select: dropped one variable (keep)
mutate: changed 2,196 values (<1%) of 'growth_dat' (0 new NA)

Plot predictions

# Trimmed plot
nd |> 
  filter(growth_dat == "Y") |> 
  ggplot(aes(yday, y = year, fill = pred, color = pred)) +
    geom_raster() +
    facet_wrap(~area, ncol = 3) +
    scale_fill_viridis(option = "magma") +
    scale_color_viridis(option = "magma") +
    labs(x = "Yearday", y = "Year", color = "Predicted SST (°C)", fill = "Predicted SST (°C)")
filter: removed 613,782 rows (62%), 378,810 rows remaining

# Full plot
nd |>
  ggplot(aes(yday, y = year, fill = pred, color = pred)) +
    geom_raster() +
    facet_wrap(~area, ncol = 3) +
    scale_fill_viridis(option = "magma") +
    scale_color_viridis(option = "magma") +
    labs(x = "Yearday", y = "Year", color = "Predicted SST (°C)", fill = "Predicted SST (°C)")

Detailed exploration of predictions

# Loop trough all areas, plot temperature as a function of yday, color by data source, facet by year

for(i in unique(nd$area)) {
  
  plotdat <- nd |> filter(area == i)
  
  print(
    ggplot(plotdat, aes(yday, pred, color = source)) + 
      scale_color_brewer(palette = "Dark2") + 
      facet_wrap(~year) + 
      geom_point(data = filter(d, area == i & year > min(plotdat$year)), size = 0.2,
                 aes(yday, temp, color = source)) + 
      geom_line(linewidth = 0.3) + 
      labs(title = paste("Area = ", i), color = "", linetype = "") + 
      guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
      theme_sleek(base_size = 8) +
      theme(legend.position = c(0.8, 0.08), 
            legend.direction = "horizontal",
            legend.spacing.y = unit(-0.3, "cm")) + 
      labs(x = "Day of the year", y = "Predicted SST (°C)")
  )
  
  ggsave(paste0(home, "/figures/supp/full_model/temp_pred_yday_area_", i, ".pdf" ), width = 17, height = 17, units = "cm")
  
}
filter: removed 901,458 rows (91%), 91,134 rows remaining
filter: removed 94,021 rows (96%), 3,449 rows remaining
filter: removed 901,458 rows (91%), 91,134 rows remaining
filter: removed 88,009 rows (90%), 9,461 rows remaining

filter: removed 901,458 rows (91%), 91,134 rows remaining
filter: removed 85,514 rows (88%), 11,956 rows remaining

filter: removed 901,458 rows (91%), 91,134 rows remaining
filter: removed 89,320 rows (92%), 8,150 rows remaining

filter: removed 901,458 rows (91%), 91,134 rows remaining
filter: removed 87,348 rows (90%), 10,122 rows remaining

filter: removed 901,458 rows (91%), 91,134 rows remaining
filter: removed 86,597 rows (89%), 10,873 rows remaining

filter: removed 901,458 rows (91%), 91,134 rows remaining
filter: removed 92,662 rows (95%), 4,808 rows remaining

filter: removed 901,458 rows (91%), 91,134 rows remaining
filter: removed 88,609 rows (91%), 8,861 rows remaining

filter: removed 901,458 rows (91%), 91,134 rows remaining
filter: removed 94,066 rows (97%), 3,404 rows remaining

filter: removed 901,458 rows (91%), 91,134 rows remaining
filter: removed 95,353 rows (98%), 2,117 rows remaining

filter: removed 947,574 rows (95%), 45,018 rows remaining
filter: removed 88,498 rows (91%), 8,972 rows remaining

filter: removed 956,358 rows (96%), 36,234 rows remaining
filter: removed 82,293 rows (84%), 15,177 rows remaining

Area-specific models

spec_preds <- list()

# Drop VN, no logger data? 
dspec <- d |> filter(!area == "VN")
filter: removed 2,129 rows (2%), 95,341 rows remaining
for(i in unique(dspec$area)) {
  
  dd <- dspec |> filter(area == i)

  if(unique(dd$area) %in% c("BS", "BT", "FB", "FM", "MU", "SI_EK", "TH")) { # RA, JM, HO, SI_HA remains
    
    mspec <- sdmTMB(temp ~ 0 + source_f + year_f + s(yday, bs = "cc"), 
                    data = dd,
                    family = student(df = 6),
                    spatial = "off",
                    spatiotemporal = "off",
                    knots = list(yday = c(0.5, 364.5)),
                    control = sdmTMBcontrol(newton_loops = 1)) 
  
  } else {
    
    mspec <- sdmTMB(temp ~ 0 + source_f + year_f + s(yday, bs = "cc"), 
                    data = dd,
                    family = student(df = 10),
                    spatial = "off",
                    spatiotemporal = "off",
                    knots = list(yday = c(0.5, 364.5)),
                    control = sdmTMBcontrol(newton_loops = 1)) 
    
  }
  
  
  sanity(mspec)

  # QQ plots
  mcmc_res_msep <- residuals(mspec, type = "mle-mcmc",
                             mcmc_samples = sdmTMBextra::predict_mle_mcmc(mspec,
                                                                          mcmc_iter = 201,
                                                                          mcmc_warmup = 200))
  
  print(ggplot(dd, aes(sample = mcmc_res_msep)) +
    stat_qq(size = 0.75, shape = 21, fill = NA) +
    stat_qq_line() +
    labs(y = "Sample Quantiles", x = "Theoretical Quantiles", title = paste("Area = ", i)) + 
    theme(aspect.ratio = 1))
  
  ggsave(paste(home, "/figures/supp/qq_temp", i, ".pdf", sep = ""), width = 11, height = 11, units = "cm")

  # Predict on new data
  nd_area <- data.frame(expand.grid(yday = seq(min(dd$yday), max(dd$yday), by = 1),
                                    area = unique(dd$area),
                                    source = unique(dd$source_f),
                                    year = unique(dd$year))) |>
    mutate(source_f = as.factor(source),
           year_f = as.factor(year)) |> 
    left_join(gdat, by = "area") |> 
    mutate(area = as.factor(area),
           growth_dat = ifelse(year >= min & year <= max, "Y", "N"))
  
  nd_area$pred <- predict(mspec, newdata = nd_area)$est

  # Save!
  spec_preds[[i]] <- nd_area
  
}
filter: removed 91,880 rows (96%), 3,461 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf

SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000444 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.44 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 2.2592 seconds (Warm-up)
Chain 1:                0.003148 seconds (Sampling)
Chain 1:                2.26235 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     83,664
           >                 ========
           > rows total       83,664
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 86,369 rows (91%), 8,972 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.00145 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 14.5 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 5.84882 seconds (Warm-up)
Chain 1:                0.017694 seconds (Sampling)
Chain 1:                5.86651 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 42 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     46,116
           >                 ========
           > rows total       46,116
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 85,868 rows (90%), 9,473 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001138 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 11.38 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 7.12418 seconds (Warm-up)
Chain 1:                0.008437 seconds (Sampling)
Chain 1:                7.13262 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     90,885
           >                 ========
           > rows total       90,885
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 83,373 rows (87%), 11,968 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.00149 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 14.9 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 7.12698 seconds (Warm-up)
Chain 1:                0.009846 seconds (Sampling)
Chain 1:                7.13682 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     91,134
           >                 ========
           > rows total       91,134
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 87,179 rows (91%), 8,162 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.00105 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 10.5 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 4.96599 seconds (Warm-up)
Chain 1:                0.0076 seconds (Sampling)
Chain 1:                4.97359 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     90,885
           >                 ========
           > rows total       90,885
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 85,207 rows (89%), 10,134 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001311 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 13.11 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 6.70454 seconds (Warm-up)
Chain 1:                0.009923 seconds (Sampling)
Chain 1:                6.71447 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     84,660
           >                 ========
           > rows total       84,660
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 84,456 rows (89%), 10,885 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001398 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 13.98 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 5.96622 seconds (Warm-up)
Chain 1:                0.04286 seconds (Sampling)
Chain 1:                6.00908 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     91,134
           >                 ========
           > rows total       91,134
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 90,521 rows (95%), 4,820 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000619 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 6.19 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 3.31345 seconds (Warm-up)
Chain 1:                0.002263 seconds (Sampling)
Chain 1:                3.31571 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     91,134
           >                 ========
           > rows total       91,134
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 86,468 rows (91%), 8,873 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.00134 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 13.4 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 5.71711 seconds (Warm-up)
Chain 1:                0.008596 seconds (Sampling)
Chain 1:                5.72571 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     91,134
           >                 ========
           > rows total       91,134
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 80,164 rows (84%), 15,177 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.001871 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 18.71 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 9.64926 seconds (Warm-up)
Chain 1:                0.013614 seconds (Sampling)
Chain 1:                9.66287 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 50 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     54,900
           >                 ========
           > rows total       54,900
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA
filter: removed 91,925 rows (96%), 3,416 rows remaining
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf

Warning in min(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf
Warning in min(x, na.rm = na.rm): no non-missing arguments to min; returning
Inf
Warning in max(x, na.rm = na.rm): no non-missing arguments to max; returning
-Inf


SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.00042 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 4.2 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:   1 / 201 [  0%]  (Warmup)
Chain 1: Iteration:  20 / 201 [  9%]  (Warmup)
Chain 1: Iteration:  40 / 201 [ 19%]  (Warmup)
Chain 1: Iteration:  60 / 201 [ 29%]  (Warmup)
Chain 1: Iteration:  80 / 201 [ 39%]  (Warmup)
Chain 1: Iteration: 100 / 201 [ 49%]  (Warmup)
Chain 1: Iteration: 120 / 201 [ 59%]  (Warmup)
Chain 1: Iteration: 140 / 201 [ 69%]  (Warmup)
Chain 1: Iteration: 160 / 201 [ 79%]  (Warmup)
Chain 1: Iteration: 180 / 201 [ 89%]  (Warmup)
Chain 1: Iteration: 200 / 201 [ 99%]  (Warmup)
Chain 1: Iteration: 201 / 201 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 1.64832 seconds (Warm-up)
Chain 1:                0.00326 seconds (Sampling)
Chain 1:                1.65158 seconds (Total)
Chain 1: 
mutate: new variable 'source_f' (factor) with 3 unique values and 0% NA
        new variable 'year_f' (factor) with 83 unique values and 0% NA
left_join: added 2 columns (min, max)
           > rows only in x        0
           > rows only in y  (    11)
           > matched rows     90,885
           >                 ========
           > rows total       90,885
mutate: converted 'area' from character to factor (0 new NA)
        new variable 'growth_dat' (character) with 2 unique values and 0% NA

nd_area <- dplyr::bind_rows(spec_preds)

# In order to have a lower temperature in before-nuclear times (without any data to inform that), we can use the nearby areas.. so FM informs BT prior to nuclear
nd_area_sub <- nd_area |> 
  mutate(keep = "N",
         keep = ifelse(area == "FM" & year <= 1980, "Y", keep), # use FM instead of BT
         keep = ifelse(area == "SI_EK" & year <= 1972, "Y", keep)) |> # use SI_EK instead of SI_HA
  filter(keep == "Y") |> # Now change the labels to BT and SI_EK...
  mutate(area = ifelse(area == "FM", "BT", "SI_HA"))
mutate: new variable 'keep' (character) with 2 unique values and 0% NA
filter: removed 825,279 rows (91%), 81,252 rows remaining
mutate: converted 'area' from factor to character (0 new NA)
# Bind rows and plot only the temperature series we will use for growth modelling
nd_area <- bind_rows(nd_area, nd_area_sub) |>
  select(-keep) |> 
  mutate(growth_dat = ifelse(area == "SI_HA" & year %in% c(1966, 1967), "Y", growth_dat)) # SI_EK and SI_HA do not have the same starting years, so we can't use allo columns from SI_EK
select: dropped one variable (keep)
mutate: changed 2,196 values (<1%) of 'growth_dat' (0 new NA)
nd_area |> 
  filter(area %in% c("FM", "BT", "SI_EK", "SI_HA")) |> 
  filter(year <= 1980 & year >= 1966) |> 
  group_by(area, year) |> 
  summarise(mean_temp = mean(pred)) |> 
  ungroup() |> 
  pivot_wider(names_from = area, values_from = mean_temp)
filter: removed 623,247 rows (63%), 364,536 rows remaining
filter: removed 298,656 rows (82%), 65,880 rows remaining
group_by: 2 grouping variables (area, year)
summarise: now 60 rows and 3 columns, one group variable remaining (area)
ungroup: no grouping variables
pivot_wider: reorganized (area, mean_temp) into (BT, FM, SI_EK, SI_HA) [was 60x3, now 15x5]
# A tibble: 15 × 5
    year    BT    FM SI_EK SI_HA
   <dbl> <dbl> <dbl> <dbl> <dbl>
 1  1966  8.15  8.15  7.97  7.97
 2  1967  8.98  8.98  8.76  8.76
 3  1968  8.69  8.69  8.37  8.37
 4  1969  8.55  8.55  8.53  8.53
 5  1970  8.30  8.30  8.04  8.04
 6  1971  8.33  8.33  8.33  8.33
 7  1972  8.83  8.83  8.90  8.90
 8  1973  9.19  9.19  8.99 11.1 
 9  1974  9.05  9.05  8.66  8.46
10  1975  7.43  7.43  9.34 14.6 
11  1976  7.09  7.09  8.00 13.7 
12  1977  5.61  5.61  7.80 13.0 
13  1978  5.59  5.59  7.57 13.1 
14  1979  5.86  5.86  7.42 11.4 
15  1980  6.96  6.96  7.73 13.4 

Plot detailed exploration of predictions

# Loop trough all areas, plot temperature as a function of yday, color by data source, facet by year

for(i in unique(nd_area$area)) {
  
  plotdat <- nd_area |> filter(area == i)
  
  print(
    ggplot(plotdat, aes(yday, pred, color = source)) + 
      scale_color_brewer(palette = "Dark2") + 
      facet_wrap(~year) + 
      geom_point(data = filter(d, area == i & year > min(plotdat$year)), size = 0.2,
                 aes(yday, temp, color = source)) + 
      geom_line(linewidth = 0.3) + 
      labs(title = paste("Area = ", i), color = "", linetype = "") + 
      guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
      theme_sleek(base_size = 8) +
      theme(legend.position = "bottom", 
            legend.direction = "horizontal",
            legend.spacing.y = unit(-0.3, "cm")) + 
      labs(x = "Day of the year", y = "Predicted SST (°C)")
  )
  
  ggsave(paste0(home, "/figures/supp/temp_pred_yday_area_", i, ".pdf" ), width = 17, height = 17, units = "cm")
  
}
filter: removed 904,119 rows (92%), 83,664 rows remaining
filter: removed 94,021 rows (96%), 3,449 rows remaining
filter: removed 896,649 rows (91%), 91,134 rows remaining
filter: removed 88,498 rows (91%), 8,972 rows remaining

filter: removed 896,898 rows (91%), 90,885 rows remaining
filter: removed 88,009 rows (90%), 9,461 rows remaining

filter: removed 896,649 rows (91%), 91,134 rows remaining
filter: removed 85,514 rows (88%), 11,956 rows remaining

filter: removed 896,898 rows (91%), 90,885 rows remaining
filter: removed 89,320 rows (92%), 8,150 rows remaining

filter: removed 903,123 rows (91%), 84,660 rows remaining
filter: removed 87,348 rows (90%), 10,122 rows remaining

filter: removed 896,649 rows (91%), 91,134 rows remaining
filter: removed 86,597 rows (89%), 10,873 rows remaining

filter: removed 896,649 rows (91%), 91,134 rows remaining
filter: removed 92,662 rows (95%), 4,808 rows remaining

filter: removed 896,649 rows (91%), 91,134 rows remaining
filter: removed 88,609 rows (91%), 8,861 rows remaining

filter: removed 896,649 rows (91%), 91,134 rows remaining
filter: removed 82,293 rows (84%), 15,177 rows remaining

filter: removed 896,898 rows (91%), 90,885 rows remaining
filter: removed 94,066 rows (97%), 3,404 rows remaining

Compare predictions from full and area-specific models

long_preds <- bind_rows(nd_area |> filter(growth_dat == "Y") |> mutate(model = "area-specific"),
                        nd |> filter(growth_dat == "Y") |> mutate(model = "full model"))
filter: removed 560,979 rows (57%), 426,804 rows remaining
mutate: new variable 'model' (character) with one unique value and 0% NA
filter: removed 613,782 rows (62%), 378,810 rows remaining
mutate: new variable 'model' (character) with one unique value and 0% NA
for(i in unique(long_preds$area)) {
  
  plotdat <- long_preds |> filter(area == i)
  
  print(
    ggplot(plotdat, aes(yday, pred, color = source, linetype = model)) + 
      scale_color_brewer(palette = "Dark2") + 
      facet_wrap(~year) + 
      geom_point(data = filter(d, area == i & growth_dat == "Y"), size = 0.2,
                 aes(yday, temp, color = source), inherit.aes = FALSE) + 
      geom_line(linewidth = 0.3) + 
      guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
      theme_sleek(base_size = 8) +
      theme(legend.position = "bottom", 
            legend.direction = "horizontal",
            legend.spacing.y = unit(-0.3, "cm")) + 
      labs(x = "Day of the year", y = "Predicted SST (°C)",
           title = paste("Area = ", i), color = "", linetype = "")
  )
  
  ggsave(paste0(home, "/figures/supp/full_model/temp_pred_comp_area_model", i, ".pdf" ), width = 17, height = 17, units = "cm")
  
}
filter: removed 767,706 rows (95%), 37,908 rows remaining
filter: removed 95,578 rows (98%), 1,892 rows remaining
filter: removed 744,126 rows (92%), 61,488 rows remaining
filter: removed 91,098 rows (93%), 6,372 rows remaining

filter: removed 700,350 rows (87%), 105,264 rows remaining
filter: removed 88,894 rows (91%), 8,576 rows remaining

filter: removed 722,166 rows (90%), 83,448 rows remaining
filter: removed 91,175 rows (94%), 6,295 rows remaining

filter: removed 728,859 rows (90%), 76,755 rows remaining
filter: removed 90,618 rows (93%), 6,852 rows remaining

filter: removed 670,062 rows (83%), 135,552 rows remaining
filter: removed 88,543 rows (91%), 8,927 rows remaining

filter: removed 761,694 rows (95%), 43,920 rows remaining
filter: removed 95,510 rows (98%), 1,960 rows remaining

filter: removed 757,302 rows (94%), 48,312 rows remaining
filter: removed 95,234 rows (98%), 2,236 rows remaining

filter: removed 700,206 rows (87%), 105,408 rows remaining
filter: removed 90,492 rows (93%), 6,978 rows remaining

filter: removed 744,126 rows (92%), 61,488 rows remaining
filter: removed 84,366 rows (87%), 13,104 rows remaining

filter: removed 772,719 rows (96%), 32,895 rows remaining
filter: removed 95,686 rows (98%), 1,784 rows remaining

filter: removed 792,438 rows (98%), 13,176 rows remaining
filter: removed 97,092 rows (>99%), 378 rows remaining

# Join wide data to plot differences easily
wide_pred <- left_join(nd_area |> 
                         dplyr::select(yday, area, year, pred, source) |> 
                         pivot_wider(names_from = source, values_from = pred) |> 
                         rename(logger_as = logger,
                                errs_as = errs,
                                fishing_as = fishing),
                       
                       nd |> 
                         dplyr::select(yday, area, year, pred, source) |> 
                         pivot_wider(names_from = source, values_from = pred) |> 
                         rename(logger_full = logger,
                                errs_full = errs,
                                fishing_full = fishing),
                       
                       by = c("yday", "area", "year"))
pivot_wider: reorganized (pred, source) into (logger, errs, fishing) [was 987783x5, now 329261x6]
rename: renamed 3 variables (logger_as, errs_as, fishing_as)
pivot_wider: reorganized (pred, source) into (logger, errs, fishing) [was 992592x5, now 330864x6]
rename: renamed 3 variables (logger_full, errs_full, fishing_full)
left_join: added 3 columns (logger_full, errs_full, fishing_full)
           > rows only in x    33,672
           > rows only in y  ( 35,275)
           > matched rows     295,589
           >                 =========
           > rows total       329,261
diff_pred <- wide_pred |> 
  mutate(logger_diff = logger_as - logger_full,
         errs_diff = errs_as - errs_full,
         fishing_diff = fishing_as - fishing_full) |> 
  select(yday, area, year, logger_diff, errs_diff, fishing_diff) |> 
  pivot_longer(c("logger_diff", "errs_diff", "fishing_diff"))
mutate: new variable 'logger_diff' (double) with 268,506 unique values and 10% NA
        new variable 'errs_diff' (double) with 268,506 unique values and 10% NA
        new variable 'fishing_diff' (double) with 268,506 unique values and 10% NA
select: dropped 6 variables (logger_as, errs_as, fishing_as, logger_full, errs_full, …)
pivot_longer: reorganized (logger_diff, errs_diff, fishing_diff) into (name, value) [was 329261x6, now 987783x5]
for(i in unique(diff_pred$area)) {
  
  plotdat <- diff_pred |> filter(area == i)
  
  print(ggplot(plotdat, aes(yday, value, color = name)) + 
      facet_wrap(~year) +  
      scale_color_brewer(palette = "Dark2") + 
      geom_hline(yintercept = 0, linetype = 2, alpha = 0.5) + 
      geom_line() +
      labs(x = "Day of the year", y = "Predicted SST (°C)",
           title = paste("Area = ", i), color = "", linetype = ""))
}
filter: removed 904,119 rows (92%), 83,664 rows remaining
filter: removed 896,649 rows (91%), 91,134 rows remaining

Warning: Removed 46116 rows containing missing values (`geom_line()`).
filter: removed 896,898 rows (91%), 90,885 rows remaining

filter: removed 896,649 rows (91%), 91,134 rows remaining

filter: removed 896,898 rows (91%), 90,885 rows remaining

filter: removed 903,123 rows (91%), 84,660 rows remaining

filter: removed 896,649 rows (91%), 91,134 rows remaining

filter: removed 896,649 rows (91%), 91,134 rows remaining

filter: removed 896,649 rows (91%), 91,134 rows remaining

filter: removed 896,649 rows (91%), 91,134 rows remaining

Warning: Removed 54900 rows containing missing values (`geom_line()`).
filter: removed 896,898 rows (91%), 90,885 rows remaining

Plot summarized data and predictions

# Summarise data
dsum <- d |> 
  group_by(year, area, source) |> 
  summarise(temp = mean(temp)) |> 
  mutate(type = "data")
group_by: 3 grouping variables (year, area, source)
summarise: now 1,601 rows and 4 columns, 2 group variables remaining (year, area)
mutate (grouped): new variable 'type' (character) with one unique value and 0% NA
# Summarise predictions from full model and area-specific model
preds <- nd |> 
  filter(growth_dat == "Y" & source == "logger") |> 
  group_by(area, year) |> 
  summarise(temp = mean(pred)) |> 
  mutate(model = "full model")
filter: removed 866,322 rows (87%), 126,270 rows remaining
group_by: 2 grouping variables (area, year)
summarise: now 345 rows and 3 columns, one group variable remaining (area)
mutate (grouped): new variable 'model' (character) with one unique value and 0% NA
preds_area <- nd_area |> 
  filter(growth_dat == "Y" & source == "logger") |> 
  group_by(area, year) |> 
  summarise(temp = mean(pred)) |> 
  mutate(model = "area model")
filter: removed 845,515 rows (86%), 142,268 rows remaining
group_by: 2 grouping variables (area, year)
summarise: now 395 rows and 3 columns, one group variable remaining (area)
mutate (grouped): new variable 'model' (character) with one unique value and 0% NA
preds_comb <- bind_rows(preds, preds_area)

ggplot(preds_comb, aes(year, temp, color = source, linetype = model)) + 
  geom_point(data = dsum, aes(year, temp, color = source), size = 0.75, alpha = 0.75, inherit.aes = FALSE) + 
  scale_color_brewer(palette = "Accent") +
  geom_line(linewidth = 0.5, color = "grey20") + 
  facet_wrap(~area) +
  theme(legend.position = "bottom")

Make final plot using the area-specific model

# Overall t_opt from 02-fit-vbge.qmd! Update when I got final model!
topt_lwr <- 6.960337
topt <- 9.468607
topt_upr <- 16.854070

# Add latitude
area <- c("BS", "BT", "FB", "FM", "HO", "JM", "MU", "RA", "SI_EK", "SI_HA", "TH")
nareas <- length(area)
lat <- c(60, 60.4, 60.3, 60.5, 63.7, 58, 59, 65.9, 57.3, 57.4, 56.1)
lon <- c(21.5, 18.1, 19.5, 18, 20.9, 16.8, 18.1, 22.3, 16.6, 16.7, 15.9)
area_attr <- data.frame(cbind(area = area, lat = lat, lon = lon)) |>
  mutate_at(c("lat","lon"), as.numeric) |> 
  arrange(desc(lat))
mutate_at: converted 'lat' from character to double (0 new NA)
           converted 'lon' from character to double (0 new NA)
ggplot(preds_area, aes(year, temp, color = temp)) + 
  facet_wrap(~factor(area, levels = area_attr$area), ncol = 3) + 
  geom_hline(yintercept = topt_lwr, linewidth = 0.3, linetype = 2, color = "tomato3", alpha = 0.3) +
  geom_hline(yintercept = topt, linewidth = 0.3, linetype = 1, color = "tomato3", alpha = 0.3) +
  geom_hline(yintercept = topt_upr, linewidth = 0.3, linetype = 2, color = "tomato3", alpha = 0.3) +
  geom_line() +
  labs(x = "Year", y = "Model-predicted annual average temperature") + 
  scale_color_viridis(option = "magma", name = "Area") +
  guides(color = "none") 

preds |> 
  group_by(area) |> 
  summarise(min = min(year),
            max = max(year)) |> 
  arrange(min)
group_by: one grouping variable (area)
summarise: now 12 rows and 3 columns, ungrouped
# A tibble: 12 × 3
   area    min   max
   <chr> <dbl> <dbl>
 1 JM     1953  2016
 2 BT     1963  1980
 3 FM     1963  2000
 4 SI_HA  1966  1972
 5 SI_EK  1968  2015
 6 FB     1969  2016
 7 MU     1981  2000
 8 HO     1982  2016
 9 BS     1985  2002
10 RA     1985  2006
11 VN     1987  1998
12 TH     2000  2014
ggsave(paste0(home, "/figures/annual_average_temperature.pdf"), width = 17, height = 17, units = "cm")
# Save prediction df
write_csv(preds_area, paste0(home, "/output/gam_predicted_temps.csv"))

Code below is not evaluated!

Growing season? This might be different for different areas…

# Find day of the year where temperature exceeds 10C by area across years
# TODO: Also for area-specific predictions?
gs_area <- nd |> 
  group_by(area, yday) |> 
  summarise(mean_pred = mean(pred)) |>
  ungroup() |> 
  filter(mean_pred > 10) |> 
  group_by(area) |> 
  summarise(gs_min = min(yday),
            gs_max = max(yday))
group_by: 2 grouping variables (area, yday)
summarise: now 4,392 rows and 3 columns, one group variable remaining (area)
ungroup: no grouping variables
filter: removed 2,706 rows (62%), 1,686 rows remaining
group_by: one grouping variable (area)
summarise: now 12 rows and 3 columns, ungrouped
nd <- left_join(nd, gs_area, by = "area")
left_join: added 2 columns (gs_min, gs_max)
           > rows only in x         0
           > rows only in y  (      0)
           > matched rows     992,592
           >                 =========
           > rows total       992,592
gs_area$mean_pred <- 10

# Plot!
nd |> 
  #filter(growth_dat == "Y") |> 
  group_by(area, yday) |> 
  summarise(mean_pred = mean(pred)) |> 
  ggplot(aes(yday, mean_pred)) +
  geom_line() +
  labs(x = "Yearday", y = "Year", color = "Predicted SST (°C)", fill = "Predicted SST (°C)") + 
  facet_wrap(~area) +
  geom_point(data = gs_area, aes(gs_min, mean_pred), inherit.aes = FALSE, color = "tomato2") + 
  geom_point(data = gs_area, aes(gs_max, mean_pred), inherit.aes = FALSE, color = "tomato2") +
  geom_hline(yintercept = 10, linetype = 2)
group_by: 2 grouping variables (area, yday)
summarise: now 4,392 rows and 3 columns, one group variable remaining (area)

Now see if there is a systematic pattern in the difference between predicted and observed logger data, which could indicate that the source effect isn’t global but area-specific.

dlog <- d |> 
  filter(source == "logger") |> 
  mutate(type = "data",
         id = paste(area, year, yday, sep = "_")) |> 
  select(id, temp) |> 
  group_by(id) |> 
  summarise(obs = mean(temp)) # sometimes we have more than 1 observation per id
filter: removed 33,486 rows (34%), 63,984 rows remaining
mutate: changed 63,984 values (100%) of 'id' (63984 fewer NA)
        new variable 'type' (character) with one unique value and 0% NA
select: dropped 18 variables (year, area, yday, month, date, …)
group_by: one grouping variable (id)
summarise: now 60,967 rows and 2 columns, ungrouped
# dlog |> 
#   group_by(id) |> 
#   summarise(n = n()) |> 
#   distinct(n)

preds_log <- nd |> 
  filter(growth_dat == "Y" & source == "logger") |> 
  mutate(type = "model",
         id = paste(area, year, yday, sep = "_")) |> 
  filter(id %in% unique(dlog$id)) |> 
  ungroup() |> 
  left_join(dlog, by = "id")
filter: removed 866,322 rows (87%), 126,270 rows remaining
mutate: new variable 'type' (character) with one unique value and 0% NA
        new variable 'id' (character) with 126,270 unique values and 0% NA
filter: removed 97,965 rows (78%), 28,305 rows remaining
ungroup: no grouping variables
left_join: added one column (obs)
           > rows only in x        0
           > rows only in y  (32,662)
           > matched rows     28,305
           >                 ========
           > rows total       28,305
p1 <- preds_log |> 
  mutate(resid = pred - obs) |> 
  ggplot(aes(as.factor(area), resid, group = as.factor(area))) +
  #geom_jitter(alpha = 0.05, color = "grey20", height = 0, width = 0.2) + 
  geom_violin(fill = "grey70", color = NA) +
  geom_boxplot(width = 0.2, outlier.colour = NA, outlier.color = NA, outlier.fill = NA) +
  guides(color = "none") +
  geom_hline(yintercept = 0, linetype = 2, color = "tomato3", linewidth = 0.75) + 
  labs(x = "Area", y = "Manual residuals")
mutate: new variable 'resid' (double) with 28,305 unique values and 0% NA
p1

p1 + facet_wrap(~year) 

Some extra plots for BT especially, because it seems that the offset for logger isn’t that big in years without any logger data.. First predict BT only but with all data sources. Note this code is not run anymore and is instead expanded above to work on all areas

# Predict
ndbt <- data.frame(expand.grid(yday = seq(min(d$yday), max(d$yday), by = 1),
                               year = seq(1980, #filter(gdat, area == "BT")$min, 
                                          2004), #filter(gdat, area == "BT")$min),
                               source = unique(d$source))) |>
  mutate(area = "BT") |> 
  mutate(area = as.factor(area),
         source_f = as.factor(source),
         year_f = as.factor(year)) 

ndbt$pred <- predict(m, newdata = ndbt)$est

# Plot
ggplot(ndbt, aes(yday, pred, color = source)) + 
  scale_color_brewer(palette = "Accent") + 
  facet_wrap(~year) + 
  geom_point(data = filter(d, area == "BT" & year >= min(ndbt$year) & year <= max(ndbt$year)), size = 0.2,
             aes(yday, temp)) + 
  geom_line() + 
  labs(title = "Area = BT", color = "") + 
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
  theme_sleek(base_size = 8) +
  theme(legend.position = "bottom", 
        legend.direction = "horizontal",
        legend.spacing.y = unit(-0.3, "cm")) + 
  labs(x = "Day of the year", y = "Predicted SST (°C)")

ggsave(paste0(home, "/figures/supp/BT_test_plots.pdf" ), width = 17, height = 17, units = "cm")

# Right, so there is an offset... but it's not big. And the cold temperatures in
# certain years in the 1980's is not due to the offset not working, just that the year
# effect in that year is only informed by cold temperatures and there's no "memory"... 
# Though it seems also that there could be a bigger offset.. perhaps this indicates
# the source should indeed vary by area? Maybe, we if we fit models by area separately instead? Then we don't
# need the area interaction, and we can instead use a random walk or AR1 structure one
# the year effect?

dbt <- d |> filter(area == "BT")

mbt <- sdmTMB(temp ~ 0 + source_f + year_f + s(yday, bs = "cc"), 
              data = dbt,
              family = student(df = 5),
              spatial = "off",
              spatiotemporal = "off",
              knots = list(yday = c(0.5, 364.5)),
              control = sdmTMBcontrol(newton_loops = 1))

sanity(mbt)
summary(mbt)

# Predict
ndbt2 <- data.frame(expand.grid(yday = seq(min(d$yday), max(d$yday), by = 1),
                                year = seq(1980, #filter(gdat, area == "BT")$min, 
                                           2004), #filter(gdat, area == "BT")$min),
                                source = unique(d$source))) |>
  mutate(source_f = as.factor(source),
         year_f = as.factor(year)) 

ndbt2$pred <- predict(mbt, newdata = ndbt2)$est

# Plot without data
ggplot(filter(d, area == "BT" & year >= min(ndbt$year) & year <= max(ndbt$year)),
       aes(yday, temp, color = source)) + 
  geom_line(data = ndbt, aes(yday, pred, color = source, linetype = "main model"), alpha = 0.2) +
  geom_line(data = ndbt2, aes(yday, pred, color = source, linetype = "BT specific model"), alpha = 0.2) + 
  scale_color_brewer(palette = "Dark2") + 
  facet_wrap(~year) + 
  labs(title = "Area = BT", color = "", linetype = "") + 
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
  theme_sleek(base_size = 8) +
  theme(legend.position = "bottom", 
        legend.direction = "horizontal",
        legend.spacing.y = unit(-0.3, "cm")) + 
  labs(x = "Day of the year", y = "Predicted SST (°C)")

ggsave(paste0(home, "/figures/supp/BT_test_plots2.pdf" ), width = 17, height = 17, units = "cm")

# Plot with data
ggplot(filter(d, area == "BT" & year >= min(ndbt$year) & year <= max(ndbt$year)),
       aes(yday, temp, color = source)) + 
  geom_point(size = 0.001, alpha = 0.3) + 
  geom_line(data = ndbt, aes(yday, pred, color = source, linetype = "main model"), alpha = 0.7) +
  geom_line(data = ndbt2, aes(yday, pred, color = source, linetype = "BT specific model"), alpha = 0.7) + 
  scale_color_brewer(palette = "Dark2") + 
  facet_wrap(~year) + 
  labs(title = "Area = BT", color = "", linetype = "") + 
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
  theme_sleek(base_size = 8) +
  theme(legend.position = "bottom", 
        legend.direction = "horizontal",
        legend.spacing.y = unit(-0.3, "cm")) + 
  labs(x = "Day of the year", y = "Predicted SST (°C)")

ggsave(paste0(home, "/figures/supp/BT_test_plots3.pdf" ), width = 17, height = 17, units = "cm")

# Plot only new predictions, color by year facet by source?
p1 <- ggplot() + 
  geom_line(data = ndbt2, aes(yday, pred, color = year, group = year), alpha = 0.7, linewidth = 0.2) + 
  scale_color_viridis() + 
  facet_wrap(~source) + 
  labs(title = "Area = BT", color = "", linetype = "") + 
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
  theme_sleek(base_size = 8) +
  theme(legend.position = "bottom", 
        legend.direction = "horizontal",
        legend.spacing.y = unit(-0.3, "cm")) + 
  coord_cartesian(ylim = c(0, 25)) +
  labs(x = "Day of the year", y = "Predicted SST (°C)", title = "BT specific model")

p2 <- ggplot() + 
  geom_line(data = ndbt, aes(yday, pred, color = year, group = year), alpha = 0.7, linewidth = 0.2) + 
  scale_color_viridis() + 
  facet_wrap(~source) + 
  labs(title = "Area = BT", color = "", linetype = "") + 
  guides(color = guide_legend(title.position = "top", title.hjust = 0.5)) + 
  theme_sleek(base_size = 8) +
  theme(legend.position = "bottom", 
        legend.direction = "horizontal",
        legend.spacing.y = unit(-0.3, "cm")) + 
  coord_cartesian(ylim = c(0, 25)) +
  labs(x = "Day of the year", y = "Predicted SST (°C)", title = "main model")

p1 / p2

ggsave(paste0(home, "/figures/supp/BT_test_plots4.pdf" ), width = 17, height = 17, units = "cm")

# Now plot the annual means. Focus on 1980-2004, so that we don't need to worry about
# the Forsmark predictions in years before heating

ndbt_sum <- ndbt |> 
  filter(source == "logger") |> 
  group_by(year) |> 
  summarise(mean_pred = mean(pred)) |> 
  mutate(model = "main model")

ndbt2_sum <- ndbt2 |> 
  filter(source == "logger") |> 
  group_by(year) |> 
  summarise(mean_pred = mean(pred)) |> 
  mutate(model = "BT-specific model")

bt_sums <- bind_rows(ndbt_sum, ndbt2_sum)

ggplot(bt_sums, aes(year, mean_pred, color = model)) + 
  geom_line() +
  labs(x = "Year", y = "Model-predicted annual average temperature") + 
  scale_color_brewer(palette = "Accent")

ggsave(paste0(home, "/figures/supp/BT_test_plots5.pdf" ), width = 17, height = 17, units = "cm")

If we want to get uncertainty, we can use nsim instead; this simulates from the linear predictor using the inverse precision matrix, which is a fast way to get a distribution of samples from which we can take e.g. quantiles and means. However, it’s still slow, so the code below isn’t executed yet.

nd_sim <- data.frame(expand.grid(yday = seq(min(d$yday), max(d$yday), by = 1),
                                 area = unique(d$area),
                                 year = unique(d$year))) |>
  mutate(source = "logger") |>
  mutate(id = paste(year, area, sep = "_"),
         source_f = as.factor(source),
         year_f = as.factor(year))

# Trim!
nd_sim <- left_join(nd_sim, gdat, by = "area")

nd_sim <- nd_sim |>
  mutate(growth_dat = ifelse(year > min, "Y", "N")) |>
  filter(growth_dat == "Y") |>
  filter(yday %in% c(gs_min:gs_min)) |>
  mutate(area = as.factor(area))

# Predict!
nsim <- 500
m_pred_sims <- predict(m, newdata = nd_sim, nsim = nsim)

# Join sims with prediction data
nd_sim_long <- cbind(nd_sim, as.data.frame(m_pred_sims)) |>
    pivot_longer(c((ncol(nd_sim) + 1):(nsim + ncol(nd_sim))))

# Summarize sims over growing season
sum_pred_gs <- nd_sim_long |>
    ungroup() |>
    group_by(year, area) |>
    summarise(lwr = quantile(value, prob = 0.1),
              est = quantile(value, prob = 0.5),
              upr = quantile(value, prob = 0.9)) |>
    ungroup()

# In order to have a lower temperature in before-nuclear times (without any data to inform that), we can use the nearby areas..
sum_pred_gs <- preds |>
  mutate(keep = "Y",
         keep = ifelse(area == "BT" & year < 1980, "N", keep),
         keep = ifelse(area == "SI_HA" & year < 1972, "N", keep)) |>
  filter(keep == "Y")

sum_pred_gs_sub <- preds |>
  mutate(keep = "N",
         keep = ifelse(area == "FM" & year < 1980, "Y", keep), # use FM instead of BT
         keep = ifelse(area == "SI_EK" & year < 1972, "Y", keep)) |> # use SI_EK instead of SI_HA
  filter(keep == "Y")

# Now change the labels to BT and SI_EK...
sum_pred_gs_sub <- sum_pred_gs_sub |>
  mutate(area = ifelse(area == "FM", "BT", "SI_HA"))

# Bind rows and plot only the temperature series we will use for growth modelling
sum_pred_gs <- bind_rows(sum_pred_gs, sum_pred_gs_sub) |> select(-keep, -type)

order <- sum_pred_gs |>
  group_by(area) |>
  summarise(mean_temp = mean(temp)) |>
  arrange(desc(mean_temp))